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ABSTRACT 


Extensions  have  been  made  to  the  NIST  3-D  cement  hydration  and  microstructure  de- 
velopment model  to  incorporate  pozzolanic  materials  such  as  silica  fume  and  fly  ash.  These 
additions  will  enable  the  simulation  of  the  hydration  and  microstructure  development  of 
blended  cements  as  well  as  conventional  Portland  cements.  SEM  and  image  analysis  tech- 
niques originally  developed  for  differentiating  the  cement  clinker  phases  present  in  a 2-D 
image  of  cement  particles  have  been  modified  and  apphed  to  equivalent  images  of  fly  ash 
particles.  Reactions  have  been  proposed  and  the  appropriate  volume  stoichiometries  deter- 
mined for  the  various  components  of  fly  ash.  Some  preliminary  calibrations  of  reaction  rate 
constants  and  dissolution  rates  have  been  performed  by  comparison  to  experimental  results 
for  the  adiabatic  heat  signatures  of  a variety  of  concretes  with  various  water-to-cement  ratios 
and  fly  ash  contents.  Further  research  needed  to  complete  this  cahbration  is  outlined. 

Keywords:  Building  technology,  cement  hydration,  computer  modelling,  fly  ash,  microstruc- 
ture,  pozzolan,  silica  fume,  simulation. 
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1 Introduction 


This  document  describes  the  efforts  to  date  for  incorporating  fly  ash  and  silica  fume  into 
the  NIST  3-D  Portland  cement  hydration  and  microstructure  development  model  [1].  A 
series  of  fly  ashes  (from  the  U.S.  and  France)  have  been  imaged  using  a combination  of  SEM 
backscattered  electron  and  X-ray  signals  to  distinguish  the  major  phases  present  (sihca,  an 
aluminosilicate,  calcium  chloride,  tricalcium  aluminate,  a calcium  aluminosihcate,  anhydrite, 
and  inert  material).  Computational  programs  for  distributing  the  cement  and  fly  ash  phases 
amongst  a size  distribution  of  three-dimensional  particles  have  been  developed.  The  fly  ash 
components  may  be  distributed  randomly  amongst  a subset  of  the  particles,  or  as  individual 
monophase  particles.  Reactions  have  been  established  for  these  starting  materials  and  the 
phases  of  ordinary  Portland  cement  based  on  information  available  in  the  literature  [2,  3]. 
Volume  stoichiometries  and  heats  of  reaction  have  been  postulated  for  these  reactions,  which 
then  have  been  included  in  a new  version  of  the  NIST  3-D  microstructural  model.  Some 
vahdation  against  adiabatic  heat  release  curves  has  been  performed  for  one  of  the  French 
fly  ashes.  The  basic  procedures  are  thus  in  place  for  modelling  the  hydration  behavior  of  a 
variety  of  blended  cements  (and  concretes),  but  further  validation  and  adjustment  of  some 
constants  (reactions  rates,  etc.)  is  still  needed.  A suggested  approach  for  performing  this 
vahdation  is  also  presented. 

The  new  version  of  the  hydration  model  and  the  auxihary  programs  pro\dded  in  Appendix 
B are  available  \da  anommious  ftp  from  edsel.cbt.nist.gov  (129.6.104.138).  The  user  can  log 
in  as  “anon^mious”  and  supply  their  e-mail  address  as  the  password.  The  new  programs  can 
be  found  in  the  /ftp/pub/CEMHYD3D/flyash  subdirectory.  The  programs  have  been 
developed  and  executed  on  a UNIX-based  computer  system  and  should  require  no  more  than 
16  iVIbytes  of  memory  for  execution.  In  addition,  a postscript  version  of  this  documentation 
can  be  found  in  the  /ftp/pub/CEMHYD3D/flyash/manual  subdirectory.  Complete 
documentation  for  the  3-D  cement  hydration  and  microstructure  development  program  has 
been  pubhshed  pre\dously  [1]  (also  available  in  the  /ftp/pub/CEMHYD3D/manual  sub- 
directory) and  should  be  consulted  along  with  this  update. 

2 Experimental  Procedures 

SEM  and  X-ray  analysis  pro\dde  a convenient  and  powerful  means  for  examining  the  mi- 
crostructure of  fly  ash  particles  [2,  4].  Here,  the  SEM/imaging  techniques  developed  pre- 
viously [5]  for  imaging  Portland  cements  have  been  modified  for  the  analysis  of  fly  ash 
samples.  In  addition  to  the  backscattered  electron  (BSE)  image.  X-ray  images  are  collected 
for  Al,  Ca,  Cl,  S,  and  Si.  Based  on  the  intensities  of  these  signals,  a separation  into  the 
following  phases  is  performed:  silica  (Si02)^  an  aluminosilicate  (5E4/2O5),  calcium  chloride 
(CaC/2),  anhydrite  (CaSO^)^  a calcium  aluminosilicate  {C aSi2Al20s)-,  tricalcium  aluminate 
(C a^AhOo) ^ and  inert  (other)  material  [2,  3].  The  initial  segmented  image  is  typically  pro- 
cessed with  a median  filter  [1]  to  improve  the  image  quahty  by  removing  much  of  the  noise 
within  indi\ddual  phase  regions.  A series  of  final  processed  images  for  a variety  of  fly  ashes  is 
provided  in  Appendix  A of  this  document.  These  final  images  can  be  analyzed  to  determine 
phase  volume  fractions,  phase  surface  area  fractions,  and  correlation  functions  to  be  used  in 
reconstructing  representative  three-dimensional  images  [1]. 
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In  reality,  of  course,  the  phase  composition  of  the  fly  ash  is  more  complex  than  being 
simply  composed  of  the  compounds  hsted  above.  However,  in  approximating  the  fly  ash  to 
be  composed  of  these  compounds,  we  are  operating  in  much  the  same  spirit  as  that  which 
utihzes  the  Bogue  compounds  €28^  and  C4AF)  as  a representation  of  a Portland 

cement. 


3 Three-Dimensional  Reconstruction 

Once  the  phase  compositions,  phase  surface  areas,  and  correlation  functions  have  been  deter- 
mined, the  same  computational  techniques  developed  for  three-dimensional  Portland  cement 
particles  can  be  used  to  reconstruct  three-dimensional  fly  ash  particles  following  a specific 
particle  size  distribution  [1,  6].  In  fact,  using  digitized  spheres  to  represent  the  fly  ash  par- 
ticles should  in  general  be  more  reahstic  than  in  the  case  of  cement  particles,  since  fly  ash 
particles  are  generally  much  more  spherical  in  shape  than  cement  particles.  In  addition  to 
the  computer  programs  for  performing  the  reconstruction  based  on  autocorrelation  analysis, 
two  new  computer  programs  for  distributing  the  fly  ash  phases  amongst  the  fly  ash  particles 
have  been  developed.  The  first  simply  distributes  the  phases,  in  their  appropriate  volume 
fractions,  totally  at  random  amongst  all  of  the  pixels  initially  identified  as  fly  ash.  The  sec- 
ond achieves  this  random  distribution  amongst  the  fly  ash  particles,  so  that  all  particles  are 
monophase,  but  the  user-supphed  appropriate  phase  volume  fractions  are  nonetheless  main- 
tained. This  program  has  been  developed  because  it  has  been  observed  in  SEM  images  that 
the  fly  ash  particles  are  often  monophase.  Listings  for  these  two  programs  (distfarand.c 
and  distfapart.c)  are  provided  in  Appendix  B of  this  document. 

During  the  assignment  of  phases  to  the  fly  ash  particles,  the  tricalcium  aluminate  con- 
tained in  the  fly  ash  is  assigned  a different  phase  ID  than  that  used  for  the  tricalcium 
aluminate  present  in  the  Portland  cement.  However,  within  the  microstructural  model,  the 
two  are  considered  to  be  equivalent  and  participate  in  all  of  the  same  reactions,  with  the 
same  dissolution  probabilities  and  diffusion  rates.  While  fly  ash  often  replaces  cement  on  a 
mass  basis,  for  use  in  the  model,  this  substitution  must  be  converted  to  a volume  basis.  A 
specific  gravity  of  2.2  is  generally  used  for  the  fly  ash  particles  for  performing  this  conversion. 
When  using  the  program  genpartSd.c  (listing  provided  in  Appendix  B),  the  fly  ash  and 
cement  can  follow  different  particle  size  distributions,  but  the  user  should  be  sure  to  place 
all  of  the  particles  in  order  of  largest  to  smallest,  regardless  of  their  phase  identification. 


4 Additional  Reactions 

The  needed  physical  properties  of  the  phases  foimd  in  fly  ash  and  the  relevant  hydration 
products  are  summarized  in  Table  1.  While  specific  gravities,  molecular  weights,  and  molar 
volumes  are  readily  available  in  the  literature  [3,  7],  heat  of  formation  data  have  yet  to  be 
located  for  all  of  the  phases.  Figure  1 summarizes  the  proposed  reactions  between  the  fly  ash 
and  the  cement  phases  and  their  hydration  products,  such  as  CH.  The  numbers  below  each 
reaction  indicate  the  volume  stoichiometries  (on  a pixel  basis)  which  must  be  maintained 
by  the  computer  model.  Based  on  these  volume  stoichiometries,  all  but  one  of  the  reactions 

^Conventional  cement  chemistry  notation  is  used  throughout  this  document  with  C = CaO,  S = Si02, 
A = AI2O3,  F = FesOs,  H = H2O,  and  5 = SO3. 
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in  Figure  1 are  seen  to  further  contribute  to  the  chemical  shrinkage  [6]  occuring  during 
hydration,  as  the  volume  of  the  solid  hydration  products  is  less  than  that  of  the  reactants 
(solids  and  water).  The  only  exception  to  this  is  the  conversion  of  primary  C—S—H  to 
“pozzolanic”  C—S—H. 

Table  1:  Physical  Properties  of  Proposed  Components  of  Fly  Ash 


Compound 

Name 

Compound 

Formula 

Density 

(Mg/m^) 

Molar  volume 
(cm^/mol) 

Heat  of  formation 
(kJ/mol) 

Silica 

S' 

2.2 

27. 

-907.5 

Aluminosilicate 

AS 

3.247 

49.9 

Strathngite 

C2ASHS 

1.94 

215.63 

Anhydrite 

cs 

2.61 

52.16 

-1424.6 

Gypsum 

CSH2 

2.32 

74.2 

-2022.6 

Calcium  chloride 

CaCh 

2.15 

51.62 

-795.8 

FriedePs  salt 

C^A(CaCh)Hw 

1.892 

296.66 

Calcium  Hydroxide 

CH 

2.24 

33.1 

-986.1 

Calcium  sihcate  hydrate 

Ci.jSH^ 

2.12 

108. 

-3283 

Pozzolardc  C—S—H 

Ci.^SHs,s 

1.69 

101.8 

-2299.1 

Tri  calcium  aluminate 

CsA 

3.03 

89.1 

-3587.8 

Dicalcium  aluminosihcate 

C2AS 

3.05 

90. 

Calcium  aluminodisihcate 

CAS2 

2.77 

100.62 

The  new  version  of  the  model  allows  for  the  conversion  of  primary  C—S—H  to  pozzolanic 
C—S—H.,  for  the  case  where  there  is  an  excess  of  silica  due  to  the  presence  of  fly  ash  or 
sihca  fume.  The  user  can  select  to  have  this  reaction  included  with  or  eliminated  from  those 
executed  in  the  model  through  the  use  of  a global  input  flag.  This  reaction  allows  for  the 
observed  reduction  in  the  C/S  ratio  of  C—S—H  over  time  in  systems  containing  pozzolans  [8]. 
Wdien  activated,  the  global  parameter  PCSH2CSH  controls  the  probability  of  a primary 
C—S—H  pixel  being  converted  to  pozzolanic  C—S—H  during  any  given  dissolution  cycle  of 
the  model.  The  calcium  hberated  by  this  reaction  is  converted  to  diffusing  CH  species  which 
can  form  solid  CH  or  participate  in  further  pozzolanic  reactions  with  the  available  silica  and 
aluminosilicates . 

For  the  purposes  of  incorporating  fly  ash,  four  new  diffusing  species  [1]  have  been  in- 
troduced into  the  model:  diffusing  anhydrite,  diffusing  (7aC/2,  diffusing  AS,  and  diffusing 
CAS2.  The  silica  itself  does  not  dissolve,  but  diffusing  CH  species  react  at  sihca  surfaces  to 
form  pozzolanic  C—S—H.  The  interactions  between  the  new  diffusing  species  and  the  phases 
of  the  Portland  cement  are  as  fohows: 

diffusing  CaCl2'-  when  a diffusing  CaCh  species  collides  with  either  solid  or  diffusing 
C3A,  Friedel’s  salt  is  formed.  If  it  colhdes  with  solid  C4AF,  FriedePs  salt,  Cif,  and 
F Hs  are  formed. 

diffusing  AS:  when  a diffusing  AS  species  collides  with  either  solid  or  difihsing  CH., 
strathngite  is  formed. 
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Pozzolanic  C—S—H  Formation 


I.ICH  + 5 + 2m  Ci.iS'ifs.g 
1.349  1 1.867  3.771 

Ci.iSH^  + 0.5F  Ci.i5i73.9  + 0.6Ci7 

1 0.083  0.943  0.184 


Calcium  Chloride- Aluminate  Reactions 

CsA  ^CaCh^  lOF  -^C3A[CaCl2)H^o 
1.7261  1 3.487  5.747 

C^AF  + CaCh  + 14F  C3A(CaC/2)i7io  A CH  + FH3 
2.4797  1 4.882  5.747  0.641  1.3522 


Anhydrite  to  Gypsum  Conversion 

CSa2H-^  CSH2 
1 0.69  1.423 


Stratlingite  Formation 

2CHa  as  + 6H  -^C2ASHs 
1 0.7538  1.631  3.26 

CA52  + C3A  + m 2C2ASH8 

1 0.886  2.862  4.286 

CA52  + C4AF  + 20^ 2C2A5i78  + CH  A FHs 
0.786  1 2.81  3.37  0.2586  0.5453 

Figure  1:  Cement-fly  ash  model  reactions  - numbers  below  reactions  indicate  volume  stoi- 
chiometries. 
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diffusing  CAS2:  when  a diffusing  CAS2  species  collides  with  solid  or  diffusing  C3A, 
stratlingite  is  formed.  K it  collides  with  sohd  C4AF,  stratlingite,  CH^  and  F H3  are 
formed. 

diffusing  anhydrite:  when  a diffusing  anhydrite  species  collides  with  either  solid  or  dif- 
fusing gypsum,  it  is  converted  into  sohd  gypsum. 

Globally,  when  the  model  is  executed  under  adiabatic  conditions,  the  activation  energy 
for  the  fly  ash  reactions  (considered  to  be  on  the  order  of  80  kJ/mole),  is  used  to  update 
both  the  probabihty  of  a diffusing  CH  species  reacting  at  a sihca  surface  and  the  probabihty 
of  dissolution  of  the  aluminosilicate  and  calcium  aluminosihcate  components  of  the  fly  ash. 
The  reactivity  of  the  fly  ash  should  perhaps  also  be  a function  of  pH,  but  this  has  yet  to 
be  incorporated  into  the  model.  Perhaps  using  an  exponential  function  to  describe  the  pH 
evolution  over  time  or  as  a function  of  degree  of  hydration,  such  as: 

pH  = pH^nit  + 7 X [1  - exp(,  ---)]  (1) 

i — Q 

would  pro\dde  a basis  for  estimating  pH  during  the  course  of  the  hydration.  In  the  above 
equation,  pHinit  is  the  initial  pH  of  the  pore  solution  when  the  cement  and  fly  ash  are  mixed 
with  water,  7 is  the  increase  in  pH  during  the  entire  course  of  hydration,  and  a is  the  degree 
of  hydration.  The  reactivity  of  the  fly  ash  could  then  be  adjusted  based  on  the  current  pH 
of  the  pore  solution. 

The  new  version  of  the  model  can  be  found  in  the  programs  disrealnew.c  and  hydreal- 
new.c  in  the  anonymous  ftp  directory  as  described  in  section  1.  An  example  annotated 
input  dataflle  for  the  execution  of  disrealnew  is  given  in  Appendix  C. 


5 Comparison  to  Experimental  Adiabatic  Heat  Sig- 
nature Data 

For  one  of  the  French  fly  ashes,  an  extensive  experimental  data  set  of  adiabatic  heat  sig- 
natures is  available.  For  a variety  of  water- to- cement  (w/c)  ratios  by  mass,  and  fly  ash 
contents,  concretes  have  been  prepared  and  characterized  by  monitoring  their  temperature 
rise  over  time  when  hydrating  under  adiabatic  conditions.  This  data  set  can  be  used  for  a 
direct  comparison  with  model  predictions.  This  fly  ash  is  a relatively  simple  one,  in  that 
it  can  be  approximated  as  containing  only  three  components:  silica  (S),  an  aluminosihcate 
(AS),  and  inert  material.  In  addition,  from  SEM  and  image  analysis,  it  appears  that  for 
this  fly  ash,  most  of  the  particles  (unlike  typical  cement  particles)  are  monophase.  Thus,  we 
can  use  the  program  distfapart.c  in  Appendix  B to  distribute  the  fly  ash  phases  randomly 
amongst  monophase  particles  to  achieve  the  appropriate  volume  fractions  (determined  to  be 
0.33  for  AS,  0.39  for  S,  and  0.28  for  inert  material  based  on  the  oxide  composition  deter- 
mined for  the  fly  ash).  Computationally,  this  approach  is  much  easier  than  using  the  detailed 
autocorrelation  analysis  to  distribute  the  phases,  as  is  typicahy  performed  for  a cement  [1]. 

The  dissolution  probabihty  of  the  aluminosihcate  and  the  reaction  probabihty  of  the  sihca 
(parameter  PPOZZ)  were  adjusted  to  provide  a ‘Test”  fit  to  the  experimental  adiabatic  heat 
signature  data  at  w/c=0.65  and  50%  fly  ash  content.  These  values  were  then  held  constant  in 
simulations  for  other  w/c  ratios  and  fly  ash  content  combinations.  A heat  of  hydration  value 
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of  800  J/g  AS  reacted  was  used  in  the  model  for  the  aluminosilicate  to  stratlingite  conversion. 
This  value  was  determined  based  on  experiments  in  which  the  fly  ash  alone  was  hydrated  in 
a solution  of  lime  (CH)  and  the  system  temperature  rise  monitored.  An  activation  energy  of 
83.14  kJ/mole  (the  same  as  that  determined  for  the  pozzolanic  reaction  of  silica  fume  [9,  10]) 
was  used  in  the  model  for  all  of  the  reactions  involving  one  or  more  components  of  the  fly  ash. 
Example  comparisons  of  the  experimental  and  model  temperature  rise  vs.  time  curves  are 
pro\dded  in  Figures  2,  3,  and  4.  In  general,  the  agreement  between  model  and  experiment 
is  always  within  a few  degrees  Celsius,  as  had  been  observed  previously  for  systems  with 
and  without  sihca  fume  additions  [10].  The  higher  model  temperatures  at  longer  times  are 
perhaps  due  to  more  water  being  incorporated  into  the  cement-fly  ash  hydration  products 
than  that  indicated  by  the  stoichiometries  in  Figure  1.  Since  the  hydration  is  being  executed 
under  sealed  curing  conditions,  it  wiU  effectively  terminate  when  all  remaining  porosity  is 
empty  (as  opposed  to  water-filled)  due  to  chemical  shrinkage.  Another  possibility  is  that 
certain  hydration  products  (e.g.,  ettringite)  are  unstable  at  higher  temperatures.  If  their 
deterioration  reactions  are  endothermic  in  nature,  the  predicted  temperature  rise  based  on 
the  computer  model  would  exceed  that  observed  experimentally.  Because  this  fly  ash  is 
rather  simple  in  structure  and  contains  only  a few  distinct  phases,  further  validation  of  the 
3-D  microstructural  model  with  other  fly  ashes  will  definitely  still  be  needed,  as  outhned  in 
the  next  section  of  this  report. 


Figure  2:  Comparison  of  experimental  (data  points)  and  cahbrated  simulated  (solid  fines) 
adiabatic  heat  signature  curves  for  w/c=0.65  concrete  with  50%  fly  ash  content. 


6 Needed  Research 

The  best  validation  of  the  3-D  microstructural  model  incorporating  fly  ash  would  be  a com- 
parison of  the  evolution  of  the  model  volumetric  phase  fractions  over  time  with  the  equivalent 
experimental  data.  For  crystalline  phases,  quantitative  X-ray  diffraction  [11]  would  appear 
to  be  the  most  straightforward  means  for  obtaining  phase  quantification.  This  would  be 
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Time  (h) 

Figure  3:  Comparison  of  experimental  (data  points)  and  simulated  (solid  lines)  adiabatic 
heat  signature  curves  for  w/c=0.45  concrete  with  50%  fly  ash  content. 


particularly  useful  for  the  fly  ash  hydration  products  such  as  Freidel’s  salt  and  stratlingite. 
Alternately,  SEM/image  analysis  may  be  useful  in  analyzing  real  microstructures  at  various 
hydration  ages.  Thermogravimetric  analysis  provides  data  on  the  water  released  at  different 
temperatures,  from  which  one  can  quantitatively  estimate  the  CH  content  of  cement /fly  ash 
pastes.  For  example,  this  technique  has  been  successfully  used  to  determine  the  appropriate 
phase  and  volume  stoichiometries  for  the  reaction  of  silica  fume  with  cement  [12]. 

Ideally,  these  experimental  studies  should  be  conducted  imder  isothermal  conditions  at 
several  temperatures.  This  would  enable  the  determination  of  activation  energies  and  the 
variation  of  reaction  products  with  temperature.  These  studies  could  then  be  extended  to 
hydration  under  adiabatic  conditions  as  presented  in  the  pre\dous  section.  K the  phases 
forming  at  high  temperatures  are  signiflcantly  different  from  those  at  lower  temperatures, 
extensive  modiflcations  to  the  computer  codes  may  be  required.  For  example,  ettringite  is 
known  to  be  unstable  at  temperatures  above  70  °C  [3],  so  that  the  model  could  be  modified 
to  make  ettringite  soluble  whenever  the  system  temperature  exceeds  70  ®C,  in  addition  to 
the  current  criteria  for  solubility  based  on  the  fraction  of  the  initial  gypsum  remaining  in 
the  system. 
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Figure  4:  Comparison  of  experimental  (data  points)  and  simulated  (solid  lines)  adiabatic 
heat  signature  curves  for  w/c=0.30  concrete  with  20%  fly  ash  content. 
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9 Appendix  A-  Images  of  Fly  Ashes 


Figure  5:  Initial  BSE  two-dimensional  image  of  French  fly  ash.  Image  is  approximately  250 
fim  X 200  fim. 
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Figure  6:  Final  segmented  two-dimensional  image  of  French  Municipal  Solid  Waste  Incin- 
eration fly  ash.  Phases  are  color-coded  as:  red,  S:  blue,  AS;  green,  CAS2-  orange,  CaClj] 
aqua,  anh^Mrite;  and  white,  inert.  Image  is  approximately  500  pm  x 400  pm. 
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Figure  7:  Final  segmented  two-dimensional  image  of  Class  C fly  ash  containing  C3A.  Phases 
are  color-coded  as:  red.  S:  blue,  AS;  green,  CAS2-  orange.  C3A:  aqua,  anh^^drite:  and  white, 
inert.  Image  is  approximated  250  pm  x 200  pm. 
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Figure  8:  Final  segmented  two-dimensional  image  of  Class  F fly  ash.  Phases  are  color-coded 
as:  red,  S;  blue,  AS;  green,  CAS2]  orange,  C3A;  aqua,  anhydrite;  and  white,  inert.  Image  is 
approximately  250  fim  x 200  pm. 
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10  Appendix  B-  Computer  Code  Listings 

10.1  ranl.c-  Random  number  generation 

/*  Random  number  generator  rani  from  Computers  in  Physics  */ 
/*  Volume  6 No.  5,  1992,  522-524,  Press  and  Teukolsky  */ 

/*  To  generate  real  random  numbers  0. 0-1.0  */ 

/*  Should  be  seeded  with  a negative  integer  */ 

#define  lA  16807 
#define  IM  2147483647 
#define  IQ  127773 
#define  IR  2836 
#define  NTAB  32 
#define  EPS  (1.2E-07) 

#define  MAX(a,b)  (a>b)?a:b 
#define  MIN(a,b)  (a<b)?a:b 

double  ranl(idum) 
int  *idum; 

{ 

int  j ,k; 

static  int  iv [NTAB] , iy=0 ; 
void  nrerrorO  ; 

static  double  NDIV  = 1 .0/(1 .0+(IM-l .0) /NTAB) ; 
static  double  RNMX  = (1.0-EPS); 
static  double  AM  = (1.0/IM); 


if  ((*idum  <=  0)  II  (iy  ==  0))  { 

*idum  = MAX(-*idum, *idum) ; 

for(j=NTAB+7; j>=0; j — ) { 

k = *idum/IQ; 

*idum  = lA* (*idum-k*IQ) -IR*k; 
if(*idum  < 0)  *idum  +=  IM; 
if(j  < NTAB)  iv[j]  = *idum; 

} 

iy  = iv[0]  ; 

} 

k = *idum/IQ; 

*idujn  = IA*(*idum-k*IQ)-IR*k; 

if(*idum<0)  *idum  +=  IM; 

j = iy*NDIV; 

iy  = iv[j]  ; 

iv[j]  = *idum; 

return  MIN(AM*iy ,RNMX) ; 

} 
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#uiidef  lA 
#undef  IM 
#undef  IQ 
#undef  IR 
#undef  NTAB 
ttundef  EPS 
#undef  MAX 
#undef  MIN 


10.2 


genpartSd.c-  3-D  particle  generation 


/* 

/* 

/* 

/* 

/* 

/* 

/* 

/* 

/* 

/* 

/* 

/* 

/* 

/* 


Program  genpartSd.c  to  generate  three-dimensional  cement 
particles  in  a 3-D  box  with  periodic  boundaries. 

Particles  are  composed  of  either  cement  clinker  or  gypsum, 
follow  a user-specified  size  distribution,  and  can 
be  either  flocculated,  random,  or  dispersed. 

Programmer:  Dale  P.  Bentz 

Building  and  Fire  Research  Laboratory 
NIST 

Building  226  Room  B-350 
Gaithersburg,  MD  20899  USA 
(301)  975-5865  FAX:  301-990-6891 

E-mail:  dale.bentz@nist.gov 


*/ 

*/ 

*/ 

*/ 

*/ 

*/ 

*/ 

*/ 

*/ 

*/ 

*/ 

*/ 

*/ 

*/ 


/*  Modified  3/97  to  allow  placement  of  pozzolanic,  inert  and  fly  ash  particles*/ 
#include  <stdio.h> 

#include  <math.h> 


#define  SYSSIZE  100  /*  system  size  in  pixels  per  dimension  */ 

#define  MAXTRIES  15000  /*  maximum  number  of  random  tries  for  sphere  placement  */ 


/*  phase  identifiers  */ 
#define  POROSITY  0 
/*  Note  that  each  particle 


#define  GEM  100  /* 
#define  CEMID  1 /* 
#define  GYPID  5 /* 
#define  AGG  6 /* 
#define  POZZID  17  /* 
#define  INERTID  18  /* 
#define  FLYASH  25  /* 
#define  NPARTC  12000  /* 
#define  BURNT  14000  /* 
#define  NUMSIZES  30  /* 


must  have  a separate  ID  to  allow  for  flocculation  */ 
and  greater  */ 

phase  identifier  for  cement  */ 
phase  identifier  for  gypsum  */ 
phase  identifier  for  flat  aggregate  */ 
phase  identifier  for  pozzolanic  material  */ 
phase  identifier  for  inert  material  */ 
phase  identifier  for  generic  fly  ash  */ 
maximum  number  of  particles  allowed  in  box*/ 
this  value  must  be  at  least  100  > NPARTC  */ 
mcLximum  number  of  different  particle  sizes  */ 


/*  data  structure  for  clusters  to  be  used  in  flocculation  */ 
struct  cluster{ 

int  part id;  /*  index  for  particle  */ 

int  clustid;  /*  ID  for  cluster  to  which  this  particle  belongs  */ 
int  partphase;  /*  phase  identifier  for  this  particle  (CEMID  or  GYPID)*/ 
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int  x,y,z,r;  /*  particle  centroid  and  radius  in  pixels  */ 

struct  cluster  *nextpart;  /*  pointer  to  next  particle  in  cluster  ♦/ 

}; 

/*  3-D  particle  structure  (each  particle  has  own  ID)  stored  in  array  cement  */ 

/*  3-D  microstructure  is  stored  in  3-D  array  cemreal  */ 

static  unsigned  short  int  cement  [SYSSIZE+1]  [SYSSIZE+1]  [SYSSIZE+1] ; 

static  unsigned  short  int  cemreal  [SYSSIZE+1]  [SYSSIZE+1]  [SYSSIZE+1] ; 

int  npart , aggsize ; /*  global  number  of  particles  and  size  of  aggregate  */ 

int  *seed;  /*  random  number  seed-  global  */ 

int  dispdist;  /♦  dispersion  distance  in  pixels  ♦/ 

int  clusleft;  /*  number  of  clusters  in  system  */ 

float  probgyp;  /*  probability  of  gypsum  particle  instead  of  cement  ♦/ 
struct  cluster  *clust [NPARTC] ; /*  limit  of  NPARTC  particles/clusters  ♦/ 

/*  Random  number  generator  rani  from  Computers  in  Physics  */ 

/*  Volume  6 No.  5,  1992,  522-524,  Press  and  Teukolsky  ♦/ 

/*  To  generate  real  random  numbers  0. 0-1.0  */ 

/♦  Should  be  seeded  with  a negative  integer  */ 

#define  lA  16807 
#define  IM  2147483647 
#define  IQ  127773 
tdefine  IR  2836 
#define  NTAB  32 
#define  EPS  (1.2E-07) 

#define  MAX(a,b)  (a>b)?a:b 
#define  MIN(a,b)  (a<b)?a:b 

double  ranl(idum) 
int  *idum; 

/*  Calls:  no  routines  ♦/ 

/*  Called  by:  gsphere ,makef loc  ♦/ 

{ 

int  j,k; 

static  int  iv [NTAB] , iy=0; 
void  nrerrorO  ; 

static  double  NDIV  = 1 . 0/ (l . 0+(IM-l . 0) /NTAB) ; 
static  double  RNMX  = (1.0-EPS); 
static  double  AM  = (1.0/IM); 

if  ((*idum  <=  0)  I I (iy  ==  0))  { 

♦idum  = MAX(-*idum,*idum) ; 
for(j=NTAB+7;j>=0;  j— ) { 
k = *idum/IQ; 

♦idum  = lA* C+idum-k+IQ) -IR*k; 
if(*idum  < 0)  *idum  +=  IM; 


20 


if(j  < NTAB)  iv[j]  = *idTiia; 


} 

iy  = iv[0]  ; 

} 

k = *iduiii/IQ; 

*idimi  = IA*(*iduin-k*IQ)-IR*k; 

if(*iduin<0)  *idum  +=  IM; 

j = iy*NDIV; 

iy  = iv[j]  ; 

iv[j]  = *idiim; 

return  MIN(AM*iy ,RNMX) ; 

} 

#undef  lA 
#undef  IM 
#undef  IQ 
#undef  IR 
#undef  NTAB 
#undef  EPS 
#undef  MAX 
#undef  MIN 

/*  routine  to  add  a flat  plate  aggregate  in  the  microstructure  */ 
void  addaggO 

/*  Calls:  no  other  routines  */ 

/*  Called  by:  main  program  */ 

{ 

int  ix,iy,iz; 
int  agglo,agghi; 

/*  Be  sure  aggregate  size  is  an  even  integer  */ 
do{ 

printf ("Enter  thickness  of  aggregate  to  place  (even  integer)  \n") ; 
scanf  ("y.d"  ,&aggsize)  ; 
printf ("%d\n" , aggsize) ; 

} while  ((aggsizey,2)  !=0)  ; 

if (aggsize ! =0) { 

agglo=(SYSSIZE/2)-((aggsize-2)/2) ; 
agghi=(SYSSIZE/2)+(aggsize/2) ; 

/*  Aggregate  is  placed  in  yz  plane  */ 
f or (ix=agglo ; ix<=agghi ; ix++) { 
f or (iy=l ; iy<=SYSSIZE; iy++) { 
f or(iz=l ; iz<=SYSSIZE; iz++) { 

/*  Mark  aggregate  into  both  particle  and  microstructure  images  */ 
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cement  [ix]  [iy] [iz]=AGG; 
cemreal  [ix]  [iy] [iz]=AGG; 

} 

> 

> 

} 

> 

/*  routine  to  check  or  perform  placement  of  sphere  of  ID  phasein  */ 

/*  centered  at  location  (xin,yin,zin)  of  radius  radd  */ 

/*  wflg=l  check  for  fit  of  sphere  ♦/ 

/*  wflg=2  place  the  sphere  */ 

/*  phasein  and  phase2  are  phases  to  assign  to  cement  and  cemreal  images  resp.  ♦/ 
int  chksph (xin , yin , zin , radd , wf Ig , phasein , phase2) 
int  xin , y in , zin , radd , wf Ig , phasein , phase2 ; 

/*  Calls:  no  other  routines  */ 

/*  Called  by:  gsphere  */ 

{ 

int  nofits,xp,yp,zp,i, j ,k; 
float  dist ,xdist ,ydist ,zdist ,ftmp ; 

nofits=^0;  /*  Flag  indicating  if  placement  is  possible  */ 

/*  Check  all  pixels  within  the  digitized  sphere  volume  */ 
f or(i=xin-radd; ((i<=xin+radd)&&(nof its="0) ) ; i++){ 
xp=i; 

/*  use  periodic  boundary  conditions  for  sphere  placement  */ 

if(xp<l)  {xp+=SYSSIZE;} 

else  if (xp>SYSSIZE)  {xp“=SYSSIZE; } 

ftmp= (float) (i“xin) ; 

xdist=f tmp*f tmp ; 

f or ( j =y in-radd ; ( ( j <=yin+radd)&& (nof it s==0) ) ; j ++) { 

yp=j ; 

/♦  use  periodic  boundary  conditions  for  sphere  placement  */ 

if(yp<l)  {yp+=SYSSIZE;} 

else  if (yp>SYSSIZE)  {yp"=SYSSIZE; } 

ftmp= (float) (j-yin) ; 

ydist=f tmp*f tmp ; 

for (k=zin-radd; ( (k<=zin+radd)&&(nof its==0) ) ;k++){ 
zp=k; 

/*  use  periodic  boundary  conditions  for  sphere  placement  */ 

if(zp<l)  {zp+=SYSSIZE;} 

else  if (zp>SYSSIZE)  {zp“=SYSSIZE; } 

ftmp= (float) (k-zin) ; 

zdist=f tmp*f tmp ; 
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/*  Compute  distance  from  center  of  sphere  to  this  pixel  ♦/ 
dist=sqrt (xdist+ydist+zdist) ; 
if ( (dist-0 . 5)<=(f loat)radd){ 

/*  Perform  placement  */ 
if (wflg==2){ 

cement  [xp]  [yp]  [zp] =phasein; 
cemreal  [xp] [yp] [zp] =phase2 ; 

} 

/*  or  check  placement  */ 

else  if ( (wf lg==l)&&(cement  [xp]  [yp]  [zp]  ! =P0R0SITY) ){ 
nof its=l ; 

} 

} 

/♦  Check  for  overlap  with  aggregate  */ 
if  ( (wf lg==l)&&( (abs(xp-( (float) (SYSSIZE+l)/2 . 0) ) )<( (float) aggsize/2 . 0) ) ){ 

nof its=l ; 

} 

} 

} 

} 

/*  return  flag  indicating  if  sphere  will  fit  */ 
return (nof its) ; 

} 

/*  routine  to  place  spheres  of  various  sizes  and  phases  at  random  */ 

/*  locations  in  3-D  microstructure  */ 

/*  numgen  is  number  of  different  size  spheres  to  place  */ 

/*  numeach  holds  the  number  of  each  size  class  */ 

/*  sizeeach  holds  the  radius  of  each  size  class  */ 

/*  pheach  holds  the  phase  of  each  size  class  */ 
void  gsphere (numgen, numeach, sizeeach, pheach) 
int  numgen; 

long  int  numeach [NUMS I ZES] ; 

int  sizeeach [NUMS IZES] , pheach [NUMS IZES]  ; 

/*  Calls:  makesph,  rani  */ 

/*  Called  by:  create  ♦/ 

{ 

int  count , x , y , z , radius , ig , tries ,phnow ; 
long  int  jg; 
float  rx,ry ,rz,testgyp ; 
struct  cluster  *partnew; 

/*  Generate  spheres  of  each  size  class  in  turn  (largest  first)  */ 
f or (ig=0 ; ig<numgen; ig++){ 
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phnow=p]ieach  [ig]  ; /*  phase  for  this  class  ♦/ 

radius=sizeeach[ig] ; /♦  radius  for  this  class  */ 

/*  loop  for  each  sphere  in  this  size  class  */ 
for(jg=l ; jg<=nuineach[ig]  ;jg++){ 

tries=0 ; 

/*  Stop  after  MAXTRIES  random  tries  */ 
do{ 

tries+=l ; 

/♦  generate  a random  center  location  for  the  sphere  */ 
x=(int) ( (float) SYSSIZE*ranl(seed) )+l ; 
y=(int) ( (float) SYSSIZE*ranl (seed) )+l ; 
z= (int ) ( (float ) SYSSIZE*ranl (seed) ) +1 ; 

/*  See  if  the  sphere  will  fit  at  x,y,z  */ 

/*  Include  dispersion  distance  when  checking  ♦/ 
/♦to  insure  requested  separation  between  spheres  */ 
count=chksph(x,y,z,radius+dispdist , 1 ,npcLrt+CEM,  0)  ; 
if (tries>MAXTRIES){ 

printf  ("Could  not  place  sphere  /.d  after  ’/.d  random  attempts  \n"  ,npart , MAXTRIES ) ; 

exit (1)  ; 

} 

} while (count ! =0) ; 


/♦  place  the  sphere  at  x,y,z  ♦/ 
npart+=l ; 

if (npart>=NPARTC){ 

printf("Too  many  spheres  being  generated  \n") ; 
printf("User  needs  to  increase  value  of  NPARTC  at  top  of  C“Code\n")  ; 

exit (1) ; 


} 

/*  Allocate  space  for  new  particle  info  */ 
dust [npart] = (struct  cluster  ♦)malloc(sizeof (struct  cluster)); 
dust  [npart]  ->partid=npart ; 
dust  [npart]  ->clustid=npart ; 

/♦  Default  to  cement  placement  ♦/ 

dust  [npart]  ->partphase=CEMID ; 

dust  [npart]  ->x=x; 

dust  [n"^  rt]  “>y=y ; 

dust  [npart]  ->z-z ; 

dust  [npart]  ~>r=radius ; 

cluslef t+=l  ; 

if (phnow==l){ 

testgyp=ranl (seed) ; 

/♦  Do  not  use  dispersion  distance  when  placing  particle  */ 
if (test gyp >probgyp) { 

count=chksph(x,y,z,radius,2,npart+CEM-l jCEMID) ; 
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} 

else{ 


/*  Place  particle  as  gypsum  */ 
count=chksph(x,y ,z, radius ,2, npart+CEM-l jGYPID) ; 
/*  Correct  phase  ID  of  particle  */ 
dust  [npart]  ->partphase=GYPID  ; 

} 

} 

/*  place  as  inert  or  pozzolanic  material  ♦/ 
else{ 

count =chksph(x,y,z, radius, 2, npart +CEM-1 ,phnow) ; 

/*  Correct  phase  ID  of  particle  */ 
dust  [npart]  ->partphase=phnow ; 

} 


} 


dust  [npart]  ->nextpart=NULL; 


} 

} 

/♦  routine  to  obtain  user  input  and  create  a starting  microstructure  */ 
void  create  0 
/*  Calls:  gsphere  */ 

/*  Called  by:  main  program  */ 

{ 

int  numsize, sphrad  [NUMSIZES] , sphase [NUMSIZES] ; 
long  int  sphnum  [NUMSIZES] ,invall ; 
int  isphjinval; 


do{ 

printf  ("Enter  number  of  different  size  spheres  to  use(max.  is  ’/.d)  \n" , NUMSIZES) ; 
scanf  ("y,d"  ,&numsize) ; 
printf("*/od  \n"  , numsize)  ; 

}while  ( (numsize>NUMSIZES) I | (numsize<0)) ; 
do{ 

printf ("Enter  dispersion  factor  for  spheres  (0-2)  \n")  ; 
scanf  ("*/,d"  ,&dispdist)  ; 
printf  ("*/,d  \n"  ,dispdist)  ; 

}while  ( (dispdist<0) I I (dispdist>2) ) ; 
do{ 

printf ("Enter  probability  for  gypsum  particles  (0. 0-1.0)  \n") ; 
scanf  ("*/,f"  ,&probgyp) ; 
printf  ("*/.f  \n"  ,probgyp)  ; 

} while  ((probgyp<0 .0) I I (probgyp>l .0)) ; 

if ( (numsize>0)&&(numsize< (NUMSIZES+1) ) ){ 
printf ("Enter  number,  radius,  and  phase  for  each  class  (largest  radius  1st)  \n" 
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printf ("Phases  are  1-  Cement  and  gypsum,  17-  Pozzolanic,  18-  Inert  25-  Fly  Ash  \n") 

/♦  Obtain  input  for  each  size  class  of  spheres  ♦/ 
f or (isph=0 ; isph<numsize ; isph++){ 

printf ("Enter  number  of  spheres  of  class  ‘/d  \n",isph+l); 

scanf  ("/old"  ,&invall)  ; 

printf ("*/.ld  \n",invall); 

sphnum[isph] =invall ; 

do{ 

printf  ("Enter  radius  of  spheres  of  class  y«d  \n",isph+l); 
printf  ("  (Integer  <=y,d  please)  \n"  ,SYSSIZE/5)  ; 
scanf  ("V.d"  ,&inval)  ; 
printf ("'/od  \n",inval); 

} while  ((invaKl)  I I (inval>(SYSSIZE/5)))  ; 

sphrad [isph] =inval ; 

do{ 

printf  ("Enter  phase  of  spheres  of  class  ‘/od  \n",isph+l); 
scanf  ("y,d"  ,&inval)  ; 
printf ("y,d  \n",inval); 

} while  ( (inval ! =l)&&(inval ! =12)&&(inval ! = 13)&&(inval !=25)  ) 
sphase [isph] =inval ; 


} 

gsphere(numsize,sphnum, sphrad, sphase) ; 

} 

} 

/*  Routine  to  draw  a particle  during  flocculation  routine  */ 

/*  See  routine  chksph  for  definition  of  parameters  ♦/ 
void  drawf loc (xin , yin , zin , radd , phasein , phase2) 
int  xin , yin , zin , radd , phasein ,phase2 ; 

/*  Calls:  no  other  routines  ♦/ 

/*  Called  by:  makefloc  */ 

{ 

int  xp,yp,zp,i, j ,k; 

float  dist ,xdist ,ydist ,zdist ,ftmp; 

/*  Check  all  pixels  within  the  digitized  sphere  volume  */ 
f or (i=xin-radd; (i<=xin+radd) ;i++){ 
xp=i; 

/*  use  periodic  boundary  conditions  for  sphere  placement  */ 

if(xp<l)  {xp+=SYSSIZE;} 

else  if (xp>SYSSIZE)  {xp-=SYSSIZE; } 

ftmp= (float) (i-xin) ; 

xdist=f tmp*f tmp ; 

f or(j=yin-radd; (j<=yin+radd) ; j++){ 
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} 


yp=j ; 

/*  use  periodic  boundary  conditions  for  sphere  placement  */ 

if(yp<l)  {yp+=SYSSIZE;} 

else  if (yp>SYSSIZE)  {yp-=SYSSIZE; } 

ftmp= (float) (j-yin) ; 

ydist=f tmp*f tmp ; 

f or(k=zin-radd; (k<=zin+radd) ;k++){ 

'zp=k; 

/*  use  periodic  boundary  conditions  for  sphere  placement  */ 

if(zp<l)  {zp+=SYSSIZE;} 

else  if (zp>SYSSIZE)  {zp-=SYSSIZE; > 

ftmp= (float) (k-zin) ; 

zdist=f tmp*f tmp ; 

/*  Compute  distance  from  center  of  sphere  to  this  pixel  */ 
dist=sqrt (xdist+ydist+zdist) ; 
if ( (dist-0 . 5)<=(f loat)radd){ 

/*  Update  both  cement  and  cemreal  images  */ 
cement  [xp]  [yp]  [zp] =phasein; 
cemreal  [xp]  [yp]  [zp]=phase2; 

} 

} 

} 

} 


/*  Routine  to  check  particle  placement  during  flocculation  */ 

/*  for  particle  of  size  radd  centered  at  (xin,yin,zin)  */ 

/*  Returns  flag  indicating  if  placement  is  possible  */ 
int  chkf loc(xin, yin, zin, radd) 
int  xin, yin, zin, radd; 

/*  Calls:  no  other  routines  */ 

/*  Called  by:  makefloc  */ 

{ 

int  nof  it  s , xp , yp , zp , i , j , k ; 
float  dist ,xdist ,ydist ,zdist ,ftmp ; 

nofits=0;  /*  Flag  indicating  if  placement  is  possible  */ 

/*  Check  all  pixels  within  the  digitized  sphere  volume  */ 
f or (i=xin-radd; ( (i<=xin+radd)&&(nof its==0) ) ;i++){ 
xp=i; 

/♦  use  periodic  boundary  conditions  for  sphere  placement  */ 

if(xp<l)  {xp+=SYSSIZE;> 

else  if (xp>SYSSIZE)  {xp-=SYSSIZE; > 

ftmp=(float) (i-xin) ; 
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xdist=f tmp*f tmp ; 

f or( j=yin-radd ; ((j<=yin+radd)&&(nof its==0) ) ; j++){ 

yp=j; 

/*  use  periodic  boundary  conditions  for  sphere  placement  */ 

if(yp<l)  {yp+=SYSSIZE;} 

else  if (yp>SYSSIZE)  {yp-=SYSSIZE; } 

ftmp= (float) (j-yin) ; 

ydist=ftmp*ftmp; 

f or(k=zin-radd; ((k<=zin+radd)&&(nofits==0)) ;k++){ 
zp=k; 

/*  use  periodic  boundary  conditions  for  sphere  placement  */ 

if(zp<l)  {zp+=SYSSIZE;} 

else  if (zp>SYSSIZE)  {zp-=SYSSIZE; } 

ftmp= (float) (k-zin) ; 

zdist=f tmp*f tmp ; 

/*  Compute  distance  from  center  of  sphere  to  this  pixel  */ 
dist=sqrt (xdist+ydist+zdist) ; 
if ((dist-0 . 5)<=(f loat)radd){ 

if ((cement  [xp]  [yp]  [zp]  ! =P0R0SITY) ){ 

/*  Record  ID  of  particle  hit  */ 
nof its=cement  [xp]  [yp]  [zp] ; 

> 

} 

/*  Check  for  overlap  with  aggregate  */ 

if ( (abs (xp» ( (float) (SYSSIZE+1) /2 . 0) ) ) < ( (f loat ) aggsize/2 . 0) ) { 
nof its=AGG ; 

} 

} 

} 

} 

/*  return  flag  indicating  if  sphere  will  fit  */ 
return(nof its) ; 

} 

/*  routine  to  perform  flocculation  of  particles  */ 
void  makeflocO 

/*  Calls:  drawfloc,  chkfloc,  rani  ♦/ 

/*  Called  by:  main  program  ♦/ 

{ 

int  partdo ,numf loc; 

int  nstart ; 

int  nleft,ckall; 

int  xm,yTii,zm,moveran; 

int  xp,yp,zp,rp,clushit jValkeep; 

int  iclus , cluspart [NPARTC]  ; 
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struct  cluster  *parttinp,*paxtpoint ,*partkeep; 


nstart=iipart ; /*  Counter  of  number  of  floes  remaining  */ 
f or(iclus=l ; iclus<=npart ; iclus++){ 
cluspart [iclus] =iclus ; 

} 

do{ 

printf ("Enter  number  of  floes  desired  at  end  of  routine  (>0)  \n"); 
scanf  ("’/.d"  ,&numf  loc)  ; 
printf  ("*/,d\n",numfloc)  ; 

} while  (numf loc<=0) ; 

while (nst art >numfloc){ 
nleft=0 ; 

/*  Try  to  move  each  cluster  in  turn  */ 
for(iclus=l ; iclus<=npart ; iclus++){ 
if  (dust  [iclus]  ==NULL)  { 
nleft+=l ; 

} 

else{ 

xm=ym=zm=0 ; 

/*  Generate  a random  move  in  one  of  6 principal  directions  * 
moveran=6 . *ranl (seed) ; 
switch (moveran) { 

case  0: 

xm=l ; 
breaJk; 

case  1: 

xm=(-l) ; 
break; 

case  2: 

ym=l; 

break; 

case  3: 

ym=(-l) ; 
break; 

case  4: 

zm=l ; 
break; 

case  5: 

zm=(-l) ; 
break; 

default : 


} 


break; 


/*  First  erase  all  particles  in  cluster  */ 
partpoint=clust [iclus] ; 

while (partpoint ! =NULL) { 
xp=partpoint->x ; 
yp=partpoint->y ; 
zp=partpoint->z ; 
rp=partpoint“>r ; 
drawf  loc  (xp  ,5rp,zp,rp,0,0); 
partpoint=partpoint“>nextpaxt ; 

> 

ckall=0 ; 

/*  Now  try  to  draw  cluster  at  new  location  */ 
partpoint=clust [iclus]  ; 
while ( (partpoint ! =NULL) && (ckall==0) ) { 
xp =p  art  p 0 i nt  - >x+xiii ; 
yp =p  art  p o i nt  - > y +yTn ; 
zp=partpo int “>z+zm ; 
rp=partpo int ->r ; 
ckall=chkf loc (xp , yp , zp , rp) ; 
partpoint=partpoint->nextpart ; 

} 

if (ckall==0){ 

/*  Place  cluster  particles  at  new  location  */ 
partpoint=clust [iclus] ; 
while (partpoint ! =NULL) { 

xp=partpoint->x+xm; 
yp=partpoint->y+yin; 
zp=partpoint“>z+zm; 
rp=partpoint”>r ; 
valkeep=partpoint“>partphase ; 
partdo=partpoint”>partid; 

drawf loc (xp , yp , zp , rp , partdo+CEM~ 1 , valkeep) ; 
/*  Update  particle  location  */ 
partpoint“>x=xp ; 
partpoint->y=yp ; 
partpo int “>z=zp ; 
partpoint=partpoint“>nextpart ; 

} 

} 

else{ 

/*  A cluster  or  aggregate  was  hit  */ 

/*  Draw  particles  at  old  location  ♦/ 
partpoint=clust [iclus] ; 
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/*  partkeep  stores  pointer  to  last  particle  in  list  */ 
while (partpoint ! =NULL) { 
xp=partpoint->x ; 
yp=partpoint->y; 
zp=partpoint->z ; 
rp=partpoint->r ; 
valkeep=partpoint->partphase ; 
partdo=partpoint->partid; 

drawf loc(xp,3rp,zp,rp,partdo+CEM-l , valkeep)  ; 
partkeep=partpoint ; 
partpoint=partpoint->nextpart ; 

} 

/*  Determine  the  cluster  hit  */ 
if (ckall!=AGG){ 

clushit=cluspart [ckall-CEM+1]  ; 

/*  Move  all  of  the  particles  from  cluster  clushit  to  cluster  iclus  */ 

parttmp=clust [clushit] ; 

/*  Attach  new  cluster  to  old  one  */ 
partkeep->nextpart=parttmp ; 
while  (parttmp ! =NirLL)  { 

cluspart [parttmp->partid] =iclus ; 

/*  Relabel  all  particles  added  to  this  cluster  */ 
parttmp->clustid=iclus ; 
parttmp=parttmp->nextpart ; 

} 

/*  Disengage  the  cluster  that  was  hit  */ 
dust  [clushit]  =NULL ; 
nstart-=l ; 

> 

} 

> 


} 

printf  ("Number  left  was  y.d  but  number  of  clusters  is  ’/od  \n"  ,nleft  ,nstart)  ; 
} /*  end  of  while  loop  */ 

cluslef t=nleft ; 


/*  routine  to  assess  global  phase  fractions  present  in  3-D  system  */ 
void  measure  0 

/*  Calls:  no  other  routines  */ 

/*  Called  by:  main  program  */ 

{ 

long  int  npor ,ngyp,ncem,nagg,npozz,ninert ,nflyash; 
int  i,j,k,valph; 

/*  counters  for  the  various  phase  fractions  */ 


31 


iipor=0 ; 
ngyp=0 ; 
ncein=0 ; 
nagg=0 ; 
iiinert=0 ; 
npozz=0 ; 
nf lyash=0 ; 

/*  Check  all  pixels  in  3-D  micro structure  ♦/ 
f or (i=l ; i<=SYSSIZE; i++) { 
for(j=l;j<=SYSSIZE; j++){ 
for(k=l;k<=SYSSIZE;k++){ 

valph=cemreal  [i]  [j]  [k]  ; 
if (valph==P0R0SITY)  {npor+=l;> 
else  if (valph==CEMID){ncem+=l;} 
else  if (valph==GYPID){ngyp+=l ; } 
else  if (valph==AGG)  {nagg+=l;} 
else  if (valph==POZZID)  {npozz+=l;} 
else  if (valph==INERTID)  {ninert+^l ; } 
else  if (valph==FLYASH)  {nf lyash+“l ; } 

} 

} 

} 

/♦  Output  results  */ 

printf ("Porosity=  */,ld  \n",npor); 
printf ("Cement=  %ld  \n",ncem); 
printf ("G3rpsum=  %ld  \n",ngyp); 
printf (''Aggregate=  y,ld  \n",nagg); 
printf (''Pozzolan=  ‘/.Id  \n'',npozz); 
printf  ("Inert=  /old  \n'' ,ninert)  ; 
printf ("Fly  Ash=  Xld  \n" ,nf lyash) ; 

} 

/*  Routine  to  measure  phase  fractions  as  a function  of  distance  from  ♦/ 

/*  aggregate  surface  */ 

void  measaggO 

/*  Calls:  no  other  routines  */ 

/*  Called  by:  main  program  */ 

{ 

int  phase  [24],ptot; 

int  icnt , ixlo , ixhi , iy , iz,phid, idist ; 

printf ("Distance  Porosity  Cement  Gypsum  Pozzolan  Inert  Fly  Ash\n' 
/♦  Increase  distance  from  aggregate  in  increments  of  one  */ 
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f or (idist=l ; idist<=(SYSSIZE-aggsize) /2 ; idist++) { 

/*  Pixel  left  of  aggregate  surface  */ 
ixlo=( (SYSSIZE-aggsize+2) /2)-idist ; 

/*  Pixel  right  of  aggregate  surface  */ 
ixhi=( (SYSSIZE+aggsize) /2)+idist ; 

/*  Initialize  phase  counts  for  this  distance  ♦/ 
f or (icnt=0 ; icnt<21 ; icnt++) { 
phase [icnt]=0; 

} 

ptot=0 ; 

/*  Check  all  pixels  which  are  this  distance  from  aggregate  surface  */ 
f or (iy=l ; iy<=SYSSIZE; iy++) { 
f or (iz=l ; iz<=SYSSIZE ; iz++) { 

phid=cemreal  [ixlo]  [iy]  [iz] ; 
ptot+=l ; 
if (phid<21){ 

phase [phid] +=1 ; 

> 

phid=cemreal  [ixhi]  [iy]  [iz] ; 
ptot+=l ; 
if (phid<21){ 

phase [phid] +=1 ; 

} 

} 

} 

/*  Output  results  for  this  distance  from  surface  */ 
printf("y,d  ‘/d  Xd  ‘/.d  Xd  Xd  ‘/odXn”  , idist , phase  [0]  , phase  [CEMID]  , phase  [GYPID]  ,p 

> 

} 

/*  routine  to  assess  the  connectivity  (percolation)  of  a single  phase  ♦/ 

/*  Two  matrices  are  used  here:  one  for  the  current  burnt  locations  */ 

/*  the  other  to  store  the  newly  found  burnt  locations  */ 

void  connect  0 

/*  Calls:  no  other  routines  */ 

/*  Called  by:  main  program  */ 

{ 

long  int  inew,ntop,nthrough,ncur,nnew,ntot ; 
int  i , j ,k,nmatx [29000] ,nmaty [29000] ,nmatz [29000] ; 
int  xcn, yen, zcn, npix, xl ,yl ,zl , igood,nnewx [29000] ,nnewy [29000] ,nnewz [29000] ; 
int  jneWjicur; 
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do{ 

printf ("Enter  phase  to  analyze  0)  pores  1)  Cement  \n"); 
scajif  ("yod",&npix)  ; 
printf("*/od  \n",npix); 

} while  ( (npix ! =0)&&(npix ! =1) ) ; 

/*  counters  for  number  of  pixels  of  phase  accessible  from  top  surface  */ 
/*  and  number  which  are  part  of  a percolated  pathway  */ 
ntop=0 ; 
nthrough=0 ; 

/*  percolation  is  assessed  from  top  to  bottom  only  */ 

/*  and  burning  algorithm  is  nonperiodic  in  x and  y directions  */ 
k=l; 

f or (i=l ; i<=SYSSIZE; i++) { 

for(j=l;j<=SYSSIZE;j++){ 
ncur=0 ; 
ntot=0 ; 

igood=0;  /*  Indicates  if  bottom  has  been  reached  */ 
if (((cement  [i]  [j]  [k] ==npix)&&( (cement  [i]  [j]  [SYSSIZE] ==npix) I I 
(cement  [i]  [j]  [SYSSIZE]==(npix+BURNT)))) I I ((cement  [i]  [j]  [SYSSIZE] >=CEM)&& 
(cement  [i]  [j]  [k]  >=CEM)&& (cement  [i]  [j]  [k]  <BURNT)&&(npix==l)))-[ 

/*  Start  a burn  front  */ 
cement  [i]  [j]  [k]+=BURNT; 

ntot+=l ; 
ncur+=l ; 

/*  burn  front  is  stored  in  matrices  nmat*  */ 

/♦  and  nnew*  */ 
nmatx [ncur] =i ; 
nmaty [ncur] ; 
nmatz[ncur]=l; 

/*  Burn  as  long  as  new  (fuel)  pixels  are  found  */ 
do{ 

nnew=0 ; 

f or(inew=l ; inew<=ncur; inew++){ 
xcn=nmatx [inew] ; 
ycn=nmaty [inew] ; 
zcn=nmatz [inew] ; 

/*  Check  all  six  neighbors  */ 
for(jnew=l ; jnew<=6; jnew++){ 
xl=xcn; 
yl=ycn; 
zl=zcn; 
if (jnew==l){ 

xl-=l ; 
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if (xl<l){ 

xl+=SYSSIZE; 

} 

> 

else  if(jnew==2){ 
xl+=l; 

if (xl>SYSSIZE){ 
xl-=SYSSIZE; 

} 

> 

else  if(jnew==3){ 

yi-=i; 

if(yl<l){ 

yl+=SYSSIZE; 

} 

} 

else  if(jnew==4){ 

yi+=i; 

if (yl>SYSSIZE){ 
yl-=SYSSIZE; 

} 

} 

else  if(jnew==5){ 
zl-=l; 

} 

else  if(jnew==6){ 
zl+=l ; 

> 

/*  Nonperiodic  in  z direction  so  be  sure  to  remain  in  the  3-D  box  */ 

if ( (zl>=l)&&(zl<=SYSSIZE) ) { 

if ((cement  [xl]  [yl]  [zl]  ==npix) I I ( (cement  [xl]  [yl]  [zl]>=CEM)&& 

(cement  [xl]  [yl]  [zl]<BURNT)&&(npix==l))){ 

ntot+=l ; 

cement  [xl]  [yl]  [zl]+=BURNT; 

nnew+=l ; 

if (nnew>=29000){ 

printf ("error  in  size  of  nnew  \n") ; 

} 

nnewx [nnew] =xl ; 
nnewy [nnew] =yl ; 
nnewz [nnew] =zl ; 

/*  See  if  bottom  of  system  has  been  reached  */ 

if (zl==SYSSIZE){ 
igood=l ; 

} 
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} 

} 

} 

} 

if (nnew>0){ 

iicur=nnew ; 

/♦  update  the  burn  front  matrices  */ 

f or (icur=l ; icur<=ncur ; icur++){ 
nmatx [icur] =nnewx [icur] 
nmaty [icur] =nnewy [icur] 
nmatz [icur] =nnewz [icur] 

} 

} 

}while  (nnew>0) ; 


} 


ntop+=ntot ; 
if (igood==l){ 

nthrough+=ntot ; 

} 


printf ("Phase  ID=  ‘/.d  \n",npix); 

printf ("Number  accessible  from  top=  */,ld  \n",ntop); 

printf ("Number  contained  in  through  pathways=  %ld  \n" ,nthrough) ; 

/*  return  the  burnt  sites  to  their  original  phase  values  */ 
f or (i=l ; i<=SYSSIZE; i++) { 
for(j=l;j<=SYSSIZE; j++){ 
f or (k= 1 ; k<=SYSSIZE ; k++) { 

if (cement  [i]  [j]  [k]>=BURNT){ 

cement  [i]  [j]  [k] “=BURNT ; 

} 

} 

} 

} 

} 

/*  Routine  to  output  final  microstructure  to  file  */ 
void  outmicO 

/*  Calls:  no  other  routines  ♦/ 

/*  Called  by:  main  program  */ 

{ 


FILE  *outf ile , *partf ile ; 
char  f ilen[80] ,f ilepart [80] ; 
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int  ix,iy,iz,valout ; 


printf ("Enter  name  of  file  to  save  microstructure  to  \n") ; 
scanf  ("*/,s"  ,f  ilen)  ; 
printf  ("‘/osXn"  jfilen)  ; 

outf ile=f open(f ilen, "w") ; 

printf ("Enter  name  of  file  to  save  particle  IDs  to  \n"); 
scanf  ("*/,s",filepart)  ; 
printf  ("*/,s\n"  ,filepart)  ; 

partf ile=f open(f ilepart , "w") ; 

f or(iz=l ; iz<=SYSSIZE; iz++){ 
f or (iy=l ; iy<=SYSSIZE; iy++) { 
f or (ix=l ; ix<=SYSSIZE; ix++) { 

valout=cemreal [ix] [iy]  [iz] ; 
f printf  (outf  ile , "y,ld\n"  ,valout)  ; 
valout=cement [ix] [iy] [iz] ; 
if (valout<0) {valout=0 ; } 
f printf  (partf  ile,  "y,d\n"  ,valout)  ; 

} 

} 

} 

f close(outf ile) ; 
f close(partf ile)  ; 

} 

mainO  { 

int  userc;  /*  User  choice  from  menu  */ 

int  nseed,ig, jg,kg; 

printf ("Enter  random  number  seed  value  \n") ; 
scanf  ("y.d"  ,&nseed) ; 
printf ("yd  \n",nseed); 
seed=(&nseed) ; 

/*  Initialize  counters  and  system  parameters  */ 
npart=0 ; 
aggsize=0; 
clusleft=0 ; 

/*  clear  the  3-D  system  to  all  porosity  to  start  */ 
f or (ig=l ; ig<=SYSSIZE ; ig++) { 
for(jg=l; jg<=SYSSIZE; jg++){ 


37 


for (kg= 1 ; kg<=SYSSIZE ; kg++)  { 

cement  [ig]  [jg]  [kg] =PDRGSITY; 
cemreal  [ig]  [jg]  [kg] =P0R0SITY; 

} 

} 

} 

/*  present  menu  and  execute  user  choice  */ 
do{ 

printfC”  \n  Input  User  Choice  \n'')  ; 
printf(''l)  Exit  \n'')  ; 

printf("2)  Add  spherical  particles  (cement , gypsum,  pozzolans,  etc.)  to  microstructu 

printf(''3)  Flocculate  system  by  reducing  number  of  particle  clusters  \n")  ; 

printf(''4)  Measure  phase  fractions  \n'')  ; 

printf("5)  Add  an  aggregate  to  the  microstructure  \n") ; 

printf("6)  Measure  single  phase  connectivity  (pores  or  solids)  \n'')  ; 

printf(''7)  Measure  phase  fractions  vs.  distance  from  aggregate  \n'')  ; 

printf("8)  Output  microstructure  to  file  \n'')  ; 

scanf  (''*/,d"  ,&userc) ; 
printf("%d  \n",userc); 
ff lush(stdout) ; 

switch  (userc)  { 
case  2: 

create 0 ; 
break; 

case  3: 

makeflocO  ; 
brecik; 

case  4: 

measure 0 ; 
break; 

case  5 : 

addaggO  ; 
break; 

case  6: 

connect  0 ; 
break; 

case  7 : 

if (aggsize ! =0) { 
measaggO  ; 

} 

else{ 

printfC'No  aggregate  present.  \n")  ; 

} 
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break; 


case  8: 

outmicO 

break; 

default : 


} 

} while  (userc!=l); 


break; 


10.3  distfarand.c-  Random  pixel  distribution  of  fly  ash 


#iiiclude  <stdio.h> 

#include  <math.h> 

/♦###################################################################*/ 
/*  Program  distfarand  */ 

/*  Purpose:  To  distribute  fly  ash  phases  randomly  amongst  particles  */ 
/*  Programmer:  Dale  P.  Bentz  */ 

/♦  */ 

/♦###################################################################*/ 

/♦  Define  phase  IDs  */ 

#define  FLYASH  25 
#define  C3A  35 
#define  CACL2  12 
#define  ASGID  14 
#define  CAS2ID  16 
#define  INERTID  18 
#define  POZZID  17 
#define  ANHYDRITE  20 
#define  SYSSIZE  100 

int  *seed; 

#include  "ranl.c" 

mainO 

{ 

FILE  *inf ile,*outf ile; 

char  f ilein [80] ,f ileout  [80] ; 

int  ix,iy,iz,valin,valout ; 

float  probasg ,probcacl2 ,probsio2 ; 

float  prph , probanh , probcas2 , probc3a; 

int  nseed; 

printf ("Enter  random  number  seed  value  (<0)\n"); 
scanf  ("’/od"  ,&nseed)  ; 
printf  ("yod\n"  , nseed)  ; 
seed=(&nseed) ; 

printf ("Enter  name  of  file  for  input  \n"); 
scanf  ("yoS"  ,f  ilein)  ; 
printf ("ysXn" ,f ilein) ; 
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printf ("Enter  name  of  file  for  output  \n") ; 
scanf  ("*/oS"  ,f  ileout)  ; 
printf  ("*/.s\n" ,fileout)  ; 

/*  Get  user  input  for  phase  probabilities  (volume  fractions)  */ 
printf ("Enter  probability  for  fly  ash  to  be  aluminosilicate  glass  \n") ; 
scanf  ("®/,f"  ,&probasg) ; 
printf  ("*/of\n"  ,probasg)  ; 

printf ("Enter  probability  for  fly  ash  to  be  calcium  aluminodisilicate  \n"); 
scanf  ("*/,f",&probcas2) ; 
printf  ("*/,f\n" ,probcas2)  ; 

printf ("Enter  probability  for  fly  ash  to  be  tricalcium  aluminate  \n")  ; 
scanf  ("*/of"  ,&probc3a)  ; 
printf  ("*/,f\n"  ,probc3a)  ; 

printf ("Enter  probability  for  fly  ash  to  be  calcium  chloride  \n"); 
scanf  ("*/,f"  ,&probcacl2)  ; 
printf  ("*/«f\n"  ,probcacl2)  ; 

printf  ("Enter  probability  for  fly  ash  to  be  silica  \n")  ; 
scanf  ("*/,f" ,&probsio2)  ; 
printf  ("*/,f\n" ,probsio2)  ; 

printf ("Enter  probability  for  fly  ash  to  be  anhydrite  \n"); 
scanf  ("*/,f"  ,&probanh)  ; 
printf  ("*/,f\n"  ,probanh)  ; 

inf ile=f open (filein, "r") ; 
outfile=fopen(f ileout , "w") ; 

/*  Convert  to  cumulative  probabilties  ♦/ 

/*  Order  must  match  that  used  in  for  loops  below  ♦/ 

probcacl2+=probasg ; 

probsio2+=probcacl2 ; 

probcas2+=probsio2 ; 

probanh+=probcas2 ; 

probc3a+=probanh ; 

/*  Scan  each  pixel  and  convert  all  fly  ash  pixels  to  some  */ 

/*  component  phase  */ 
f or (ix=0 ; ix<SYSSIZE; ix++) { 
f or (iy=0 ; iy<SYSSIZE ; iy++) { 
f or (iz=0 ; iz<SYSSIZE; iz++) { 

f scanf  (inf  ile,  "*/,d"  ,&valin)  ; 

valout=valin; 

if (valin==FLYASH) { 

valout=INERTID ; 
prph=ranl (seed) ; 
if (prph<probasg) { 
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valout=ASGID ; 


} 


} 

else  if (prph<probcacl2){ 
valout=CACL2 ; 

} 

else  if (prph<probsio2){ 
valout=POZZID ; 

} 

else  if (prph<probcas2){ 
valout=CAS2ID; 

} 

else  if  (prph<probaiih)  { 

val out = ANHYDRITE; 

} 

else  if (prph<probc3a){ 
valout=C3A; 

} 


fprintf  (outf  ile , ''*/,d\n"  , valout)  ; 


f close(inf ile) ; 
f close(outf ile) ; 
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10.4  distfapart.c-  Random  particle  distribution  of  fly  ash 


#iiiclude  <stdio.h> 

#include  <mat]i.h> 

/*###################################################################*/ 

/♦  Program  distf apart  to  distribute  fly  ash  phases  randomly  amongst  ♦/ 


/*  monophase  particles  */ 
/*  Programmer:  Dale  P.  Bentz  */ 
/*  Date:  May  1997  */ 
/*  */ 


/♦###################################################################*/ 

/*  Define  phase  identifiers  for  fly  ash  components  */ 

#define  FLYASH  25 
#define  C3A  35 
#define  CACL2  12 
#define  ASGID  14 
#define  INERTID  18 
#define  POZZID  17 
#define  ANHYDRITE  20 
#define  CAS2ID  16 
#define  SYSSIZE  100 
#define  NPARTC  12000 

int  *seed; 

#include  "ranl.c" 

mainO 

{ 

FILE  *inf ile,*partfile,*outf ile; 
char  f ilein [80] ,f ileout [80] ,f ilepart [80]  ; 
int  ix, iy , iZjValin, valout ,partin; 
float  probasg,probcacl2 ,probsio2; 
float  probc3a,prph,probcas2,probanh; 
static  int  phase [NPARTC] , part id  [NPARTC] ; 
int  nseed, count , cl ,phnew; 

long  int  totcnt , ascnt , cacl2cnt ,pozzcnt , inertcnt ; 
long  int  anhcnt ,cas2cnt ,c3acnt ; 

long  int  markc3a,markas ,markcacl2 ,markpozz ,markinert ,markanh,markcas2 ; 

ascnt=cacl2cnt=pozzcnt=inertcnt=0 ; 
cas2cnt=anhcnt=c3acnt=0 ; 

printf ("Enter  random  number  seed  value  (<0)\n"); 
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scanf  ("*/ocl'' ,&nseed)  ; 
printf  ("yod\ii"  ,nseed) ; 
seed=(&nseed) ; 


printf ("Enter  name  of  file  for  input  \n") ; 
scanf  ("y«s"  ,f  ilein)  ; 
printf  ("‘/osXn"  jfilein)  ; 

printf ("Enter  name  of  file  for  particle  IDs  \n"); 
scanf  ("y,s"  ,f  ilepart)  ; 
printf ("ysXn" ,f ilepart) ; 

printf ("Enter  name  of  file  for  output  \n") ; 
scanf  ("y.s"  ,f  ileout)  ; 
printf ("XsXn" ,f ileout) ; 


printf ("Enter  total  number  of  fly  ash  pixels  \n") ; 
scanf  ("y.ld"  ,&totcnt)  ; 
printf ("%ld\n" , totcnt) ; 


/*  Get  user  input  for  phase  probabilities  (volume  fractions  */ 
printf ("Enter  probability  for  fly  ash  to  be  aluminosilicate  glass  \n") ; 
scanf ("Xf" ,&probasg) ; 
printf  ("y,f\n"  ,probasg)  ; 

printf ("Enter  probability  for  fly  ash  to  be  calcium  aluminodisilicate  \n"); 
scanf  ("yof"  ,&probcas2)  ; 
printf  ("y,f\n"  ,probcas2)  ; 

printf ("Enter  probability  for  fly  ash  to  be  tricalcium  aluminate  \n")  ; 
scanf  ("yof"  ,&probc3a)  ; 
printf ("%f\n" ,probc3a) ; 

printf ("Enter  probability  for  fly  ash  to  be  calcium  chloride  \n") ; 
scanf ("%f" ,&probcacl2) ; 
printf  ("yof\n"  ,probcacl2) ; 

printf  ("Enter  probability  for  fly  ash  to  be  silica  \n"); 
scanf ("yf" ,&probsio2) ; 
printf ("%f\n" ,probsio2) ; 

printf ("Enter  probability  for  fly  ash  to  be  anhydrite  \n"); 
scanf  ("y.f"  ,&probanh)  ; 
printf  ("‘/of  \n"  ,probanh)  ; 


/*  Determine  goal  counts  for  each  phase  */ 
markas=(long) (probasg* (float) totcnt) ; 
markpozz=(long) (probsio2* (float )totcnt) ; 
markcacl2=(long) (probcacl2*(float) totcnt) ; 
mark anh= (long) (probanh* (float) totcnt) ; 
markcas2=(long) (probcas2* (float )totcnt) ; 
markc3a=(long) (probc3a* (float) totcnt) ; 
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niarkinert=(long)  ( (1 . -probasg-probsio2-probcacl2-probaiih- 
probcas2-probc3a)* (float) totcnt) ; 

/♦  Convert  probabilities  to  cumulative  */ 

/*  Order  must  be  the  scime  as  in  for  loop  below  */ 

probcacl2+=probasg ; 

probsio2+=probcacl2 ; 

probcLnh+=probsio2 ; 

probcas2+=probanh ; 

probc3a+=probcas2 ; 

inf ile=f open(f ilein, "r") ; 
partf ile=f open(f ilepaxt , "r") ; 

f or(ix=0 ; ix<NPARTC;  ix++)-[ 

phase [ix] =part id [ix] =0 ; 

} 

count=0 ; 

/*  First  scan — find  each  particle  and  assign  phases  ♦/ 
f or (ix=0 ; ix<SYSSIZE;  ix++)-C 
f or (iy=0 ; iy<SYSSIZE ; iy++) { 
f or (iz=0 ; iz<SYSSIZE; iz++) { 

f scanf  (partf  ile/'/od"  ,&partin)  ; 
f scanf  (inf  ile,  ''*/,d"  ,&valin)  ; 
if ( (valin==FLYASH)&&(partid[partin] ==0) ){ 
count+=l ; 

part id [part in] =count ; 

valout=INERTID ; 
do{ 

prph=rcLnl  (seed)  ; 

if ( (prph<probasg)&&(ascnt<markas) ) { 
valout=ASGID ; 

} 

else  if ( (prph<probcacl2)&& (cacl2cnt<markcacl2) ) { 
valout=CACL2 ; 

} 

else  if ( (prph<probsio2)&&(pozzcnt<markpozz) ) { 
valout=P0ZZID ; 

} 

else  if ( (prph<probanh)&&(anhcnt<markanh) ){ 
valout=ANHYDRITE; 

} 

else  if ( (prph<probcas2)&&(cas2cnt<markcas2) ){ 
valout=CAS2ID; 
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> 

else  if ( (prph<probc3a)&&(c3acnt<markc3a) ){ 
valout=C3A; 

> 

}wh.ile ( (valout==INERTID)&& (inertciit>markiiiert ) ) ; 

phase [count] =valout ; 

} 

if (valin==FLYASH){ 

cl=partid [partin] ; 

phnew=phase [cl] ; 

if (phnew==ASGID){ 
ascnt+=l ; 

} 

else  if (phnew==CACL2){ 
cacl2cnt+=l ; 

} 

else  if (phnew==POZZID) { 
pozzcnt+=l ; 

} 

else  if (phnew==ANHYDRITE){ 
anhcnt+=l ; 

} 

else  if (phnew==CAS2ID){ 
cas2cnt+=l ; 

} 

else  if (phnew==C3A){ 
c3acnt+=l ; 

} 

else  if (phnew==INERTID){ 


} 

> 

} 

f close(inf ile) ; 
f closeCpartf ile) ; 

/*  Now  distribute  phases  in  second  scan  ♦/ 

inf ile=f open (filein, "r”) ; 

partf ile=f open(f ilepaxt , "r") ; 

outf ile=f open(f ileout , "w”)  ; 

f or (ix=0 ; ix<SYSSIZE; ix++) { 

f or (iy=0 ; iy<SYSSIZE; iy++) { 

f or (iz=0 ; iz<SYSSIZE; iz++) { 
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f scanf  (partf  ile , ’'7od"  ,&paxtin) ; 
f scanf  (inf  ile,  "yod"  ,&valin) ; 
valout=valin; 
if (valin==FLYASH){ 

count=partid [partin] ; 
valout=phase [count] ; 

} 

fprintf  (outf  ile , "/.dXn"  , valout) ; 


} 

> 

} 

f close(inf ile) ; 
f closeCpartf ile) ; 
f close(outf ile) ; 
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new.c 

-283 

0 

cimwf aOlg . img 
1 3 4 2 5 6 
35 

pcimwf aOl . img 
0 

5000 

1 

500 

0.01  9000. 
0.0002  100000. 
0.2  2500. 

200 

5002 

6.0 

25. 

45.7 

83.14 

0.00127 

0.000 

0 

0 


ix  C-  Example  input  datafile  for  disreal- 


random  number  seed 

flag  indicating  that  final  microstructure  is  not  to  be  saved 

filename  for  file  containing  input  3-D  microstructure 

phase  IDs  assigned  to  C3S,  C2S,  C3A,  C4AF,  gypsum,  and  agg. 

phase  ID  assigned  to  C3A  in  fly  ash  particles 

filename  for  file  containing  particle  ID  microstructure 

number  of  one  pixel  particles  to  add 

number  of  cycles  of  hydration  model  to  execute 

flag  indicating  hydration  is  to  be  under  sealed  conditions 

maximum  number  of  diffuion  steps  per  cycle 

parameters  for  CH  nucleation 

parameters  for  C3AH6  nucleation 

parameters  for  FH3  nucleation 

frequency  in  cycles  to  check  pore  space  percolation 

frequency  in  cycles  to  check  total  solids  percolation 

induction  time  in  hours 

initial  temperature  in  degrees  Celsius 

activation  energy  for  cement  hydration  in  kJ/mole 

activation  energy  for  pozzolanic  reactions  in  kJ/mole 

conversion  factor  to  go  from  cycles  to  time 

mass  fraction  of  aggregates  in  mixture  proportions 

flag  indicating  hydration  is  under  isothermal  conditions 

flag  indicating  primary  to  pozzolanic  CSH  conversion  is  off 


new.c 

-283 

0 

cimwf aOlg . img 
1 3 4 2 5 6 
35 

pcimwf aOl . img 
0 

5000 

1 

500 

0.01  9000. 
0.0002  100000. 
0.2  2500. 

200 

5002 

6.0 

25. 

45.7 

83.14 

0.00127 

0.000 

0 

0 
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